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Abstract. The stellar population synthesis in unresolved composite objects is a very tricky problem. Indeed, it is a degenerate 
problem since many parameters affect the observables. The stellar population synthesis issue thus deserves a deep and rigourous 
analysis. In this paper we present a method of inversion which uses as observables the intensities at each pixel of a galactic 
spectrum and provides the stellar contribution to luminosity of all stars considered in a database. The main contribution of this 
paper to the synthesis problem is that it provides an analytical computation of the uncertainties accompanying a solution. This 
constitutes an important improvement relative to previous methods which do not provide such infomation except in the method 
described by Pelat (1997) and Moultaka & Pelat (2000). The latter uses the equivalent widths and intensities of stellar spectra 
in order to reproduce the equivalent widths of a galactic spectrum. The novelty of this work relative to the previous one is that 
the dust emission present in the IR spectra can be modeled as well as the velocity dispersion of stars that broadens the lines of 
a galactic spectrum. Tests are also performed in order to estimate the reliability of the method and the influence on the results 
of an additive continuum present in a studied spectrum, for example in the case of AGNs. 
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1. Introduction 

Thanks to the analysis made by Scheiner (1899) and later by 
Ohman (1934) who noticed the similarity between the spectra 
of the sun and the nucleus of M32, the idea of considering a 
galactic spectrum as the product of a combination of stellar 
spectra was born. Since that time, different works were carried 
out with the aim of providing the stellar populations of 
composite objects by means of stellar spectra. 
The main problem of the stellar population synthesis principle 
is that it is a degenerate problem. A large number of free 
parameters contribute to the making of the observed composite 
object spectrum (a galaxy or a cluster). Examples of such 
parameters are the stellar luminosity classes, stellar spectral 
types, metallicities, the physical characteristics and amount of 
dust, the geometry of the system, the presence of an AGN etc... 
All these parameters differently affect the observed spectrum 
of a composite object according to the considered wavelength 
domain of the spectrum (e.g. the production of absorption 
lines, emission lines, additive featureless continuum, redden- 
ing etc.). 

In order to deduce the stellar population of a galaxy or a 
cluster, two approaches have been adopted, the "direct" 
approach and the "inverse" one. In the direct approach (e.g. 
Tinsley 1972, Chariot & Bruzual 1991, Bruzual & Chariot 



1993, Leitherer et al. 1999, Fioc & Rocca-Volmerange 1997, 
Vazdekis, 1999, Bruzual & Chariot 2003), an a priori model of 
the studied object is constructed by choosing the shape of the 
initial mass function, the star formation rate, the geometry, the 
rate of supernovae and/or any other parameters. As the system 
evolves with time, synthetic observables are constructed at 
each time step with the use of stellar evolutionary models 
(the "Padova" and "Geneva" groups, e.g. Girardi et al. 2004, 
Charbonnel et al. 1999 and references herein) and directly 
compared to the real observables at all time steps. The retained 
solution is the one that best reproduces the observables. 
This approach makes it very hard to estimate the uncertainties 
on the stellar population and all the free parameters. 
In the inverse approach (e.g. Faber 1972, Joly 1974, O'Connell 
1976, Bica 1988, Schmidt et al. 1989, Silva 1991, Pelat 
1997,1998, Moultaka & Pelat 2000) one does not assume any 
predefined model for the stellar composition of the studied ob- 
ject. The observables and a method of inversion are only used 
in order to derive the characteristics of the stellar population. 
This is done by minimising a function called the objective 
function; it represents the goodness of the fit between the 
synthetic model and the observables. The estimation of the 
error bars on the stellar population is more confortable in this 
approach. However, it is usually made using Monte-Carlo 



2 



J. Moultaka: A new inverse method for stellar population synthesis and error analysis 



simulations which require a large number of simulating 
spectra. The synthesis procedure becomes thus quickly very 
costly in computing time and provides biased results. 

An inverse method described in Pelat (1997), Moultaka & 
Pelat (2000) and Moultaka et al. (2004) deals with the problem 
of minimisation in the sense that it provides a unique global 
minimum (Pelat 1997). It provides as well an analytical error 
analysis (Moultaka & Pelat 2000) and a study of the influence 
of constraints added a priori while searching for a solution 
(Moultaka et al. 2004). The method uses the equivalent widths 
(EWs) of the absorption lines of a galactic (or any composite 
object) spectrum as observables to be fitted by a combination 
of EWs of the same absorption lines and continuum intensities 
of stellar spectra. 

This method is limited as it requires the measurement of equiv- 
alent widths and consequently, an assumption of the location 
of the continuum in the spectra. Moreover, at high spectral 
resolution, the equivalent widths measurements become very 
problematic as one needs to distinguish all available lines and 
determine their intervals. On the other hand, the dust emission 
occuring particularly in the IR part of the spectrum is not 
modeled, nor is the velocity dispersion of stars which broadens 
the absorption lines and becomes very effective at high spectral 
resolution. 

In this paper, we present a new method of inversion which 
uses all spectral pixel intensities of a spectrum as observables 
instead of equivalent widths of absorption lines and we provide 
an analytical error analysis. The reddening in the optical do- 
main due to the presence of dust, the dust emission in the IR 
and the velocity dispersion of stars are also modeled. 
The method and the error analysis are described in the next sec- 
tion and tested in Sect. 4. In Sect. 3, we discuss the properties 
of the database used in the synthesis problem. An application 
to the globular cluster G 170 is described in Sec. 5. In appen- 
dices B and C are listed respectively the facilities of the code 
and the output files. 

2. The inverse method 

2.1. Search for the solution : 

First, we consider a simple case where a galaxy (or a composite 
object) is only composed of stars. This assumption leads us to 
consider the spectrum of a galaxy as the result of the sum of 
stellar spectra. Each pixel intensity of the galactic spectrum is 
then obtained as follows: 
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where Fx ga i is the galactic intensity at wavelength A and 
the intensity at the same wavelength of a star of class i (i.e. 
spectral type and luminosity class) of the stellar database. 
is the number of stars considered in the database. 

Let us call the contribution to the galactic luminosity 
at a wavelength A of stars of class i (i.e. kxi = tt^-). If we 

' Asal 



make the assumption that the galaxy is only composed of stars 
considered in the database then : 



At a reference wavelength Aq one can write , 
implies that: 

Fxgal 



•1q ga' 



(2) 



which 



A a gal 



;=1 A oi 



(3) 



or equivalently I ga ij = kjji, where k\ = k^ 0l and I ga ij and 
Iji are respectively the galactic and class i stellar intensities at 
wavelength A (indexed with j) normalised to the intensity at Ao 
(i.e. /, fl ,; = ^ and = 

A synthetic spectrum can then be constructed in a similar 
way (equivalent to Eq. Q) as follows: 



(4) 



where F Asyn is the synthetic intensity at wavelegnth A. Using 
similar notations as previously, we obtain the synthesis equa- 
tion: 



Isyn j — ^ kJj, 



(5) 



In order to ensure that a solution is physically acceptable, the 
stellar contributions to luminosity ki have to be all positive. 
Thus, considering in addition the condition of Eq. 0, the syn- 
thesis problem is finally stated as follows: 



Ilk — I s yn — 1 gal 

k >0 
bk = 1 



(6) 



In the previous system, // is the n/xn+ matrix composed of the 
n\ normalised intensities at each pixel of the n+ stellar spectra 
(i.e. composed of the 7 /( ), I ga i and I syn are respectively the vec- 
tors composed of the n\ galactic and synthetic normalised inten- 
sities at each pixel (i.e. respectiveley I ga ij and I sy „j). Finally, k 
is the vector composed of the n* stellar contributions to lumi- 
nosity kj. 

Solving the synthesis problem becomes thus to find the set k of 
stellar contributions to luminosity at the reference wavelength 
that minimises the following objective function: 



/(*) = H/ga/-//*H 2 



(7) 



This expression is equivalent to the "synthetic distance" de- 
fined by Pelat (1997) for the EW case and represents the "mean 
residual intensities". We keep the same terminology and ex- 
press the synthetic distance in our case: 



D 2 =f(k) =Tj = rihaij-lZ= l hh) 2 



— Sfelftfl/j hynj) 2 



(8) 



As shown in Eq. (|8), the synthetic distance for the case of in- 
tensities happens to be similar to the "principal synthetic dis- 
tance" defined by Pelat (1997) for the case of EWs. Therefore, 
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the problem is simplified in this case and possesses a unique 
solution. It is the solution of the least square problem subject 
to the contraints shown in equation system 

The minimisation procedure is done using the constrained 
least square method of type LS 1 of the NAG library. 

2.2. Reddening, dust emission and velocity dispersion: 

Due to the presence of dust in the Milky Way and inside the 
studied galaxy, a galactic spectrum is reddened in the optical. 
This is due to the absorption by dust of the light which is higher 
in the blue than in the red part of the spectrum. 
The intensities of a reddened spectrum are transformed as fol- 
lows: 



reddened 



(A) = F(^)10 



(-0A*F Law (A)E(B-V)) 



(9) 



where F(A) is the intensity at wavelength A, F reddenediX) the 
intensity at the same wavelength of the reddened spectrum 
and Fig W {A) the reddening law. The latter can be chosen and 
included in the code very easily (see Appendix B). The present 
law used for the synthesis in Sects. 3 and 4 is Howarth's law 
(Howarth 1983). 

In the present synthesis method, the insertion of the reddening 
effect into an optical spectrum is made by constructing differ- 
ent databases of stellar spectra reddened by a certain amount 
of E(B-V) and performing the synthesis using the different 
databases. The solution providing the smallest synthetic 
distance is retained with the corresponding value of E(B-V). 

A second effect due to the presence of dust is an additive 
emission produced by the heated dust. It can be described by a 
blackbody continuum of low temperature (T~ 50QK - \QQQK) 
and therefore, is assumed to be negligible in the optical 
domain. 

Dust emission in the IR case is included as an additional star. 
The search for a solution is made by adding to the database 
a blackbody continuum with a given temperature (or any 
other chosen function, see Appendix B) and performing the 
inversion. The test is repeated for different temperatures of 
the blackbody and, as for the reddening, the retained solution 
is the one which minimises the synthetic distance. The corre- 
sponding dust temperature is adopted for the model. 

The velocity dispersion of stars in a galaxy influences the 
spectrum of the latter by a broadening effect of the stellar ab- 
sorption lines. This effect obviously influences the results es- 
pecially at high spectral resolution. 

The velocity dispersion is included in the code by testing the re- 
sults obtained with different stellar databases constructed from 
the initial database by smoothing all the spectra with a given 
value of velocity dispersion. The convolution is made with a 
Gaussian function but any other function can be inserted easily 
in the code (see Appendix B). 



2.3. Error analysis : 

In the following, we assume that the contribution to the errors 
from the spectra of the database is negligible compared to that 



of the galaxy. If this would not be the case, there would be no 
sense in performing the synthesis with such bad data. 

Given an approximate solution I sy „o, equation system is 
still valid when I ga i is replaced by I sv „o- 
After differentiating the obtained equations, one gets the fol- 
lowing system: 



Hdk — dl syn 
bdk = 



(10) 



As the synthetic distance is a Euclidean distance, one can write 
dlsyn = Hdl ga i where H is the orthogonal projector on the hy- 
perplane tangent to the "synthetic surface" (i.e. the set of in- 
tensities having an exact solution, see Appendix A and Pelat 
1997) at /,,.„(). It follows from equation system fllOt : 



Hdk = Hdl ga i 

bdk = 



This translates into: 

dk = (ii +t n + r l n +, i Hd ^ al 



(ii) 



(12) 



where // + = 



The variance-covariance matrix of the stellar contributions (V*) 
can be deduced from Eq. fl!2i and written as a function of the 
variance-covariance matrix of the galactic normalised intensi- 
ties V/ as follows : 



V k =< dk, dk' > 

= (// + 7/ + )- 1 // + 'H < dl+dt' > H7/(// + 7/ + ) 1 



(// + '//+)- 1 // +f I | // /; // . 



gal' u gal 

ii + (n +t n + y 



o o 



(13) 



In order to compute the orthogonal projector H, the follow- 
ing two methods can be used: 

- From the first equation of equation system JlOi (Ildk = 
dl syn ) it follows that the vector of the deviation around the 
synthetic normalised intensities dl syn belongs to the linear 
envelope of the columns of matrix //. This implies that the 
linearly independent column vectors of this matrix form a 
basis of the tangent hyperplane when the origin of the vec- 
tor space of the normalised intensities is shifted to I syn o 
(because, in this case, dl syn would be on the tangent plane). 
A QR decomposition of matrix // (e.g. Golub & Van Loan 
1996) allows one to build the projector H. 

- The second way of computing the orthogonal projector H 
is to search for the equations describing the synthetic sur- 
face (Qr(Jsyni>Isyn2< ■■■Jsynr) = 0) as well as the gradients 
at l sy „o of functions Q r relative to I syn j- These gradients 
are then arranged in column vectors of a matrix G. As in 
the previous method, the orthogonal projector is deduced 
from a QR decomposition of matrix G. 

The search for functions Q r and their gradients is described 
in Appendix A. 
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Using the orthogonal projector H, one can also compute the 
error on the synthetic distance in the same way as described 
in Moultaka et al. (2004) for the EW case. A similar result is 
obtained in the case of intensities : 

a-jf. = 2 ^(I gal0 - / !V „„) r (H - ld)V/(H r - ld)(I gaI0 - /. 5V „o)(14) 

The error bars on the reddening and velocity dispersion are 
estimated using the different values of these parameters tested 
in each run and providing synthetic distances which lie inside 
the error bars of the minimal distance. 

2.4. Constraining the solution 

The solution provided by the minimisation procedure is, in 
mathematical terms, the best solution. However, this solution 
may not satisfy (astro)physical constraints which are not taken 
into account in Sect. |2] (except the positivity constraint). For 
this reason, it is crucial to constrain the solution in such a way 
that it is (astro )physically acceptable. 

This is done by minimising, as previously, the synthetic dis- 
tance defined in Eq. (|8j while the contributions to luminosity 
at the reference wavelength k are subject to linear constraints 
written as follows: 

/<{ c \}<« (15) 

In the previous equation, / and u are respectively lower 
and upper constant vector limits, C is the constraint matrix 
composed of the linear coefficients of the stellar contributions 
to luminosity k,. In order to take into account the second 
constraint shown in equation system (|6jl (bk = 1), the first 
line of matrix C is composed of a series of 1 and the second 
components of vectors / and u are equal to 1. The positivity 
constraint is included in Eq. dl5> if the first components of 
vectors / and u are respectively equal to zero and infinity. 
Note that in the present case (the case of intensities), one 
does not need to iterate on the synthetic surface as is done 
in Moultaka et al. (2004) where the observables are the EWs 
of the absorption lines, because the synthetic distance here is 
exaclty the "principal synthetic distance". 
As is described in Moultaka et al. (2004), the constraints are 
chosen in such a way that the Initial Mass Function (IMF) of 
each episode of star formation in the studied galaxy satisfies 
the standard IMFs known in the litterature. 
Two types of constraints can be used in the present version 
of the method: the "Decreasing IMF mode" where the only 
hypothesis on the IMF is that it is a decreasing function. The 
second type of constraint is the "Standard mode" where the 
constraints on the IMF are generalised in such a way that 
the standard IMFs (Salpeter 1955, Scalo 1986, Kroupa et al. 
1993) are all satisfied. The constraints are then imposed on the 
numbers (or equivalently on the contributions to luminosity) 
of dwarf stars on the one hand and on dwarf and more evolved 
stars (supergiants included) on the other hand. A detailed 
description of the constraint model is given in Moultaka et al. 
(2004). 



The facilities of the code and the description of the output 
files ar listed in Appendices B and C. 

3. Properties of the database 

The choice of the stellar (or cluster) database is of great impor- 
tance. Indeed, the solution of a synthesis depends closely on 
this choice and consequently, the interpretation of the stellar 
population of a composite object may be biased. Therefore, a 
database should be the most informative available. This means 
that the wavelength domain should be the most extended, the 
spectral resolution the highest, and the number of stars (or clus- 
ters) the largest possible, including all stellar classes (i.e. spec- 
tral types, luminosity classes and metallicities) in the case of 
stars or a large sample of ages and metallicities in the case of 
clusters. However, the spectral range and spectral resolution are 
finite and must be comparable to those of the composite object 
spectrum, which restricts the information extracted from the 
spectra. Consequently, the different components of the database 
are less distinguished and their number can hardly be as large 
as the variety of known stellar (or cluster) classes. For this rea- 
son, at a given spectral range and resolution, a database should 
be built in such a way that its different components are not lin- 
early correlated with each other or equivalently in our case, that 
no spectrum of the database can be synthesised by the remain- 
ing spectra. Indeed, if two spectra were not distinguished in the 
sense mentioned above and contribute to the synthesis solution 
of a composite object, it would be impossible to assert whether 
both stellar classes are present in the stellar population of the 
composite object or only one of them. Consequently, the num- 
ber of stars (or clusters) in the database should decrease with 
the spectral resolution or the spectral range. 
In addition, as it has been shown in Moultaka & Pelat (2000) 
(see their Fig. 2), the choice of the database has to be adapted 
to the quality of the composite object observation i.e. if two un- 
corrected stellar (or cluster) spectra were separated from each 
other within the noise level of the composite object spectrum 
(the S/N ratio of the stellar (or cluster) spectra is assumed to be 
always much higher than that of the composite object spectrum, 
see Sect. 12. 3> . they should not be both kept in the database be- 
cause the quality of the composite object spectrum would not 
allow us to distinguish between them. 

Moreover, it is obvious that the relative flux calibration of the 
stellar (or cluster) spectra is absolutely necessary and crucial 
in this method as the solutions depend closely on the differ- 
ent relative fluxes of each spectrum (see Eq. |5j». However, the 
absolute flux calibration is not needed as the spectra are all nor- 
malised to one at the reference wavelength (see Sect. 12. It but 
the magnitudes of the stellar classes would be necessary if one 
wishes to derive the number of stars from the stellar contribu- 
tions to luminosity (£,) and constrain the solutions as described 
in Sect. EH 

4. Test of the method 

A normal test of the method is to simulate a galaxy by choos- 
ing different contributions fc, of the stars as an input (kj„ put ) and 
construct the spectrum of a galaxy with the given contributions 
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and the stellar spectra. The inversion of the simulated galaxy 
with the method described in Sect. [2] will provide an estimate 
of the efficiency of the method for different values of the signal 
to noise ratio of the spectra. 

We have chosen 11 stars out of the 31 stars of Boisson et al. 
(2000) in order to perform the simulation. The 1 1 stars have 
been chosen in such a way that the different stellar classes (hot 
and cool stars, dwarfs, giants and supergiants, solar metallicity 
and SMR stars) with different continuum slopes are represented 
best. The spectral resolution of this database is of 11 A. The 
choice of this database is motivated by the fact that a higher 
number of stars would considerably increase the computing 
time for each simulation. Moreover, as our aim is to have a 
good representation of the different stellar classes, the spectral 
resolution of this database is sufficient to perform the simula- 
tions. 

The simulated galaxy has been constructed after convolving the 
stellar spectra with a Gaussian distribution of parameters fi = 
and cr = 140 km/s. The obtained spectrum is then reddened 
by E(B - V) = 0.2 mag using the reddening law of Howarth 
(1983). 

4.1. Simulating the errors 

Considering the Signal to Noise ratio (S /N) at each pixel as 
the ratio of the mean value of the normalised intensity F mean 
over its standard deviation cr F , the standard deviation of the 
normalised intensity at each pixel is obtained for each S/N ratio 
through the relation: 



We apply to the normalised intensities at each pixel of the 
simulated galaxy a random error following a normal law of pa- 
rameters fi equal to the simulated value F mean and cr to the stan- 
dard deviation calculated in Eq. dl6l for a given spectral reso- 
lution. 

Monte-Carlo simulations are then performed in order to deter- 
mine the mean values of the stellar contributions, the reddening 
and the velocity dispersion as well as their standard deviations. 

4.1.1. Results 

The results of the Monte-Carlo simulations are shown in Tables 
Hto|3]for three values of the Signal to Noise ratio (S/N ~ 25, 
50 and 100) and for three different populations : a population 
dominated by intermediate and late-type stars in Table ^ by 
young and early-type stars in Tableland by intermediate stars 
inTableE] 

These tables show clearly that for an infinite value of the S/N, 
the method provides the exact population. When the value of 
this ratio is larger than 50, the results are quite well-defined 
(with more than about 2cr error) in the example of the interme- 
diate population (Table QJ. They are also well-determined for 
the contributions exceeding 7% in the example of the interme- 
diate and late-type population (Tabled and for contributions 
exceeding 9% in the example of young early-type population 
(Table|2j. 



4.2. Simulating an additive continuum 

In order to evaluate the influence of an additive continuum to 
the solutions obtained with this method, we constructed a set of 
galaxies without any reddening or velocity dispersion, in which 
10% of the V-band luminosity arises from a power law contin- 
uum with index a = —1.5, -0.5,0.5, 1.5. Another set of galax- 
ies in which the contribution to luminosity of the continuum is 
of 50% is also constructed. The synthesis procedure is then per- 
formed on the spectrum of the constructed galaxies. The test is 
made on three different populations as described in Sect. 14.1.11 
(a population dominated by young and early-type stars, a sec- 
ond population dominated by intermediate-type stars and the 
last one dominated by late-type stars). We show the results for 
the intermediate population in Tables|4]and[5] 
The striking result is that in all cases the solutions show a large 
increase of the contributions to luminosity of hot (O- and A- 
type) stars. This result suggests that whatever the slope of the 
additive continuum (i.e. increasing or decreasing), the remedy 
of the absence of this continuum in the database is a high con- 
tribution of hot stars since these are the only spectral types 
showing a quasi-featureless continuum in the optical band. 
The high contribution of hot O- and A-type stars results in 
higher values of the reddening, especially in the case of pos- 
itive indices (increasing function) of the power law. This effect 
is expected, since the continuum shape of the hot stars is de- 
creasing, therefore, to straighten up the spectra, a higher red- 
dening is needed. 

It is also noteworthy that for a contribution to luminosity of 
50% of the additive continuum, high velocity dispersions are 
required as well to synthesise the spectrum (see Table[5}. 
This test shows that the stellar population composition is highly 
affected by the presence of an additive continuum (e.g. the pres- 
ence of an AGN or heated dust in the IR case) if this continuum 
is not modeled in the synthesis. In particular, a considerable ar- 
tificial contribution to the luminosity of young hot population 
appears in such cases. 

Moreover, a clear correlation appears between the reddening, 
the velocity dispersion and the presence of a featureless contin- 
uum. This correlation can be expected (but not easily proved) 
since the additive continuum changes the slope of the spec- 
trum and partly affects the lines in a similar way to the velocity 
dispersion. Indeed, an additive continuum dilutes the absorp- 
tion lines of a spectrum resulting in the reduction of the line 
strengths as is the case when the spectrum is smoothed. 

5. Application to the spectrum of the globular 
cluster G 170 

The method is applied to the optical spectrum of the globular 
cluster G 170 of the galaxy M31; the spectral resolution is 
of A/1 = 17 A. The database used for the synthesis is the one 
described in Boisson et al. (2000) with a spectral resolution 
of 1 1 A. The choice of this example is motivated by the fact 
that the stellar population of this cluster has been analysed in 
Moultaka et al. (2004) using the EW method and in Jablonka et 
al. (1992). It is thus an ideal example to study with the present 
method since we wish to compare the results with those of 
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star 


kl„pi,t(%) 


kouipi,A%) 


kompui(%) 


kout P ut(%) 


koutpul(%) 






S/N = oo 


S/N = 100 


S/N = 50 


S/N = 25 


O7-B0V 


1.51 


1.51 


1.46 ± 1.52 


1.55 ± 2.11 


1.67 ± 3.09 


A1-3V 


3.03 


3.03 


2.57 ± 2.33 


1.85 ±2.74 


1.87 ±3.58 


F2V 


4.54 


4.54 


4.97 ±5.12 


5.90 ± 7.04 


6.19 ± 8.56 


rG5IV 


6.06 


6.06 


6.25 ± 4.00 


6.58 ± 6.06 


6.04 ± 7.85 


rK3V 


7.58 


7.58 


7.31 ± 1.97 


6.72 ±3.56 


6.36 ±5.13 


rMlV 


9.09 


9.09 


9.07 ±0.51 


9.08 ± 1.02 


9.23 ± 2.09 


G9III 


10.61 


10.61 


11.08 ±5.85 


12.07 ± 9.99 


14.26 ± 15.89 


rK3III 


12.12 


12.12 


11.90 ±3.52 


12.09 ±6.51 


13.59 ± 10.07 


M5III 


13.64 


13.63 


13.64 ±0.19 


13.63 ±0.36 


13.68 ±0.71 


rG2I 


15.15 


15.15 


15.15 ± 3.87 


14.12 ±7.38 


11.26 ± 11.35 


M2I 


16.67 


16.67 


16.59 ±0.79 


16.41 ± 1.55 


15.86 ±3.54 


D 2 




210~ 15 


0.175 ±0.009 


0.702 ± 0.040 


2.812 ±0.172 


E(B-V) 


0.2 


0.2 


0.20 ± 0.005 


0.20 ± 0.01 


0.20 ± 0.02 


VDispersioAkm/ S) 


140 


140 


140.2 ± 6.0 


140.6 ± 11.9 


137.0 ±26.0 



Table 1. Synthesis of a simulated galaxy constructed with the contributions to luminosity listed in the second column. The stellar 
population is dominated by intermediate- and late-type stars. The different solutions correspond to different values of the S/N 
ratio (oo, 100, 50 and 25). D 1 is the synthetic distance, E(B-V) the reddening and cr the dispersion. 



star 


klnput(%) 


koutput(%) 


koutput(%) 


kou,pul(%) 


koutput(%) 






S/N = oo 


S/N= 100 


S/N = 50 


S/N = 25 


O7-B0V 


16.67 


16.67 


16.38 ± 2.09 


15.76 ±4.19 


13.34 ± 7.57 


A1-3V 


15.15 


15.15 


15.69 ±4.05 


15.55 ±6.34 


15.74 ±9.48 


F2V 


13.64 


13.64 


12.66 ± 7.47 


12.75 ± 10.82 


12.14 ± 13.14 


rG5IV 


12.12 


12.12 


12.16 ±4.02 


11.68 ±6.98 


9.94 ±9.41 


rK3V 


10.61 


10.61 


10.39 ± 1.54 


9.79 ± 3.08 


8.52 ± 5.22 


rMlV 


9.09 


9.09 


9.07 ± 0.35 


9.12 ±0.73 


9.31 ± 1.53 


G9III 


7.58 


7.58 


8.58 ± 6.05 


9.95 ± 9.43 


13.39 ± 14.55 


rK3III 


6.06 


6.06 


5.75 ±3.42 


5.97 ± 5.30 


5.77 ±6.61 


M5III 


4.54 


4.54 


4.56 ±0.11 


4.59 ± 0.21 


4.71 ±0.51 


rG2I 


3.03 


3.03 


3.29 ±2.56 


3.48 ±4.16 


5.47 ± 7.73 


M2I 


1.51 


1.51 


1.45 ±0.52 


1.35 ±0.95 


1.66 ± 1.92 


D 2 




5 10- 16 


0.008 ± 0.004 


0.340 ±0.019 


1.359 ±0.076 


E(B-V) 


0.2 


0.2 


0.20 ± 0.01 


0.20 ± 0.02 


0.19 ±0.05 


& Dispermanikm/ s) 


140 


140 


141.4 ±9.1 


139.0 ±21.9 


139.2 ±42.9 



Table 2. Same as in Tabled for a population dominated by young stars. 



star 


ki npu ,{%) 


koutput(%) 


kout P i,t(%) 


koutpiii(%) 


koutput(%) 






S/N = oo 


S/N = 100 


S/N = 50 


S/N = 25 


O7-B0V 


6.06 


6.06 


5.76 ±2.08 


5.42 ± 3.26 


5.31 ±4.81 


A1-3V 


7.58 


7.58 


7.86 ±3.80 


7.63 ± 5.30 


8.73 ± 7.42 


F2V 


10.61 


10.61 


9.95 ± 7.00 


10.47 ± 10.49 


10.36+ 9.91 


rG5IV 


16.67 


16.67 


16.83 ± 3.46 


15.90 ±6.47 


14.54 ±9.93 


rK3V 


15.15 


15.15 


15.00 ± 1.53 


14.66 ± 3.03 


14.85 ±4.52 


rMlV 


1.51 


1.51 


1.44 ±0.35 


1.48 ±0.73 


1.63 ± 1.49 


G9III 


13.64 


13.64 


14.59 ±6.28 


15.81 ± 11.31 


14.69 ± 11.93 


rK3III 


12.12 


12.12 


11.91 ± 3.61 


11.84 ± 6.51 


12.82 ±7.82 


M5III 


3.03 


3.03 


2.69 ±0.10 


2.72 ± 0.20 


3.06 ± 0.42 


rG2I 


9.09 


9.09 


9.51 ±3.11 


9.69 ± 5.86 


9.67 ± 9.22 


M2I 


4.54 


4.55 


4.45 ±0.51 


4.37pml.04 


4.33 ± 1.91 


D 2 




4 10 16 


0.086 ± 0.004 


0.344 ±0.018 


1.387 ±0.074 


E(B-V) 


0.2 


0.2 


0.20 ± 0.01 


0.20 ± 0.02 


0.20 ± 0.03 


0~Dis P ersion(km/ s) 


140 


140 


140.8 ± 6.9 


142.0 ± 13.2 


144.0 ± 20.3 



Table 3. Same as in Tabled for a population dominated by intermediate-type stars. 
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k,npuA%) 


kout P i,t(%) 


kout P ut(%) 


kout P ut(%) 


koutpi<t(%) 






a = -1.5 


-0.5 


0.5 


1.5 


hot stars 


12 


19 


21 


28 


38 


population 


inter 


inter 


inter 


inter 


inter+ young 


Continuum 


10 










D 1 




210~ 3 


3 10- 3 


6 lO" 3 


2 10~ 2 


E(B-V) 


0.0 


0.03 


0.07 


0.13 


0.23 


& Dispersion (kill/ s) 


















Table 4. Solutions of a simulated galaxy dominated by an intermediate population with a contribution to luminosity of 12% from 
hot stars and a 10% contribution from a featureless continuum (kj nput ). The different solutions correspond to the synthesis of the 
constructed galaxies where the additive continuum has a shape of a power law of different indices a. 





k[„ P ut(%) 


koutpm(%) 


koutput(%) 


koulput(%) 


koutput(%) 






a = -1.5 


-0.5 


0.5 


1.5 


hot stars 


1 


46 


47 


36 


28 


population 


inter 


young 


young 


inter+ young 


inter+young 


Continuum 


50 










D 1 




3 10 2 


410~ 2 


0.22 


1.24 


E(B-V) 


0.0 


0.17 


> 0.3 


> 0.3 


>0.3 


& Dispersion (/CfH/ S^) 





> 200 


> 200 


> 200 


> 200 



Table 5. Same as in Table|4] but the contribution to luminosity of the additive continuum is of 50%. 



the EW method and of an independent analysis. For this 
reason, the choice of the stellar database is restricted to a low 
resolution database of 31 uncorrected stars (at this spectral 
resolution and spectral range from ~ 5000A to ~ 9000A, only 
this number of stars is necessary to represent the HR diagram, 
see Sect. [3}. Obviously, the present method can be applied to 
high resolution spectra of composite objects using stars from 
the new libraries such as the Indo-US library (Valdes et al. 
2004), the UVES library (Bagnulo et al. 2003), the ELODIE 
library (Soubiran et al. 1998) or the ELODIE archive after a 
flux calibration of the spectra (Moultaka et al. 2004) and many 
others, or synthetic libraries as that of Bertone et al. (2004) 
etc... In that case, because of the higher spectral resolution, the 
computing time would increase since for the same wavelength 
range, a higher number of pixels would be present in the 
spectra and a larger number of stars would be used for the 
synthesis (because the number of uncorrected stars would also 
increase). 

The results of the synthesis obtained with the present method 
(the unconstrained and the constrained solutions with the 
"Decreasing IMF" and the "Standard" modes) are shown in 
Tables[6]and[7]and in Figs.nand[3] The results listed in Table|6] 
are obtained using all the intensities of the spectrum except 
four regions which were not corrected for the atmospheric 
bands (see Boisson et al. 2000) and the region of the Na 
absorption line (5874 to 5914A) which is also affected by the 
ISM. The spectrum of the unconstrained solution is shown in 
Fig. ^ superimposed on the dereddened observed spectrum. 
This figure shows that the spectrum is very well synthesised. 
However, the stellar population of this solution is not satisfac- 
tory because most of the contributions are not well determined 
suggesting that the stellar database is not adequate for the 



analysis of the globular cluster in that spectral range, spectral 
resolution and S/N ratio of the globular cluster spectrum. In 
addition, the contribution of hot A-type stars and the high red- 
dening of ~ 0.3 mag are probably not real since this globular 
cluster is known to be 2 10 9 yrs old (Jablonka et al. 1992) and 
reddened by ~ 0.05 mag. For this reason, we searched for 
other solutions where different parts of the spectrum were not 
used in the synthesis procedure. We found that the contribution 
of hot stars disappears when using only the red part of the 
spectrum (for instance from 6200 A to 9000 A). The results 
are shown in Table The obtained unconstrained solution is 
well defined and provides a similar stellar population as the 
one obtained by Moultaka et al. (2004). It agrees with the 
known age of this cluster but the population here is metal rich. 
The resulting value of the velocity dispersion is due to the 
non comparable spectral resolution of the stellar database in 
comparison with the globular cluster spectral resolution. The 
synthetic and observed spectra are shown in Fig. [3] 
The need of hot quasi-featureless stars in the previous solution 
obtained using the whole spectrum is thus due to the blue part 
of the spectrum. Both unconstrained solutions have similar 
populations but the solution of the synthesis of the whole 
spectrum has in addition a contribution of A-type stars because 
a bluer continuum is needed in this case. This behaviour is 
well shown in Fig. |2] where the variations of the synthetic 
distance and the contribution from hot stars with reddening 
are plotted. For all dispersions, the contribution of hot stars 
increases with the reddenning while the synthetic distance 
decreases. The optimal solution is on a border of the domain of 
parameters (at E(B-V)=0.3). This suggests that the minimum 
is a local one but a global minimum is probably at a higher 
value of the reddening which is not acceptable in physical 
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terms (for instance, the minimal solution is at E(B-V)=0.52 
mag with a contribution of 22% of hot stars). On the other 
hand, a similar plot corresponding to the synthesis of the red 
part of the spectrum (see Fig. @J shows that in this case, the 
optimal solution lies inside the domain of parameters and is 
therefore more reliable. In Moultaka et al. (2004), the solution 
does not show a contribution from hot stars. This results 
from the use of an empirical continuum to deduce the values 
of the EWs of absorption lines. The additional "empirical" 
information of the continuum used in the EWs method cannot 
be used in our case. A wider wavelength range is needed to 
provide the information about the blue part of the spectrum. 
This information is provided indirectly and empirically in 
the EWs method but is not justified as no information on the 
shape of the continuum in the blue part of the spectrum is 
given. The lack of information in the blue part of the spectrum 
is a weakness whatever the synthesis method is (direct or 
inverse). However, since in the present method the quality of 
the solution is also estimated thanks to the error analysis, it is 
possible to know the confidence level of the solution (in the 
present case for example, the solution was not well determined 
suggesting thus a lack of information). 

This analysis shows once more that the choice of the database 
and the wavelength range are crucial in the synthesis proce- 
dure and that the wavelength range should harbour enough 
signatures of the different stellar classes. It should also allow 
one to distinguish between the stellar classes within the noise 
level of the composite object spectrum as explained in Sect. [3] 

For both syntheses, the constrained "Decreasing" mode so- 
lutions are equal or equivalent to the unconstrained solution as 
was the case for the method with EWs (Moultaka et al. 2004). 
The "Decreasing IMF" mode provides solutions which are ac- 
ceptable in the sense that their synthetic distances lie inside the 
lcr error around the synthetic distance of the unconstrained so- 
lution. In Table0 the solution provides an earlier turnoff of the 
stellar population (in G4V suggesting an age of 2 10 9 yrs). 
The results are very satisfactory since the unconstrained solu- 
tions are globally the same as the constrained ones. This shows 
that the method is very stable according to such constraints i.e. 
that the unconstrained solutions are physically acceptable. 

6. Summary and conclusion 

In this paper, an inverse method of stellar population synthesis 
is described. This method uses all intensities at each pixel of 
a composite object spectrum as observables to be fitted by a 
combination of spectra from a database. 
The novelty of this work is that the reddening, velocity disper- 
sion and dust emission in the IR part of the spectrum are also 
modeled and an analytical computation of the errors around a 
solution is provided. 

Linear (astro)physical constraints can be included while 
searching for a solution. Constraints on the IMF have shown 
that the unconstrained solutions are very satisfactory as they 
satisfy such constraints. 

Tests of the method have shown that for a S/N ratio larger than 
50, the solution is well determined and when this ratio is infi- 




Fig. 2. Synthesis of the whole spectrum of the globular clus- 
ter G 170: Variation with reddening of the synthetic distance 
(top) and the contribution to luminosity of O- and A-type 
stars (k) (down) for different fixed dispersion values (<x = 
0, 20, 40, 60, 80 and lOOkm/s). The horizontal fines delimit the 
minimal synthetic distance with its lcr error. 




0.3 0.35 



Fig. 4. Same plot as in Fig.|2]but corresponding to the synthesis 
of the red part of the spectrum of G 170 (from 6200 to 9000 A). 
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G 170: Unconstrained solution 




5000 



Fig. 1. Synthetic and observed spectra respectively in dark and light lines for the unconstrained solution of the globular cluster 
G 170. 



1.2 



I 0.8 



0.6 



0.4 - 



G 170: Unconstrained solution 
~i 1 1 1 1 1 1 1 1 1 1 r 




_i i i_ 



J_ 



_i i i i_ 



J_ 



_i i i_ 



7000 



8000 



X (A) 



Fig. 3. Same plot as in Fig.^but for the synthesis performed only with the red part of the spectrum (from 6200 to 9000A). 



nite, the solution is the exact one. 

Finally, in the present spectral optical band used in this paper, it 
has been shown that the synthesis problem is very sensitive in 
the blue part. The presence of hot stars in the solutions is often 
artifactual due to a lack of information in the blue part of the 
spectrum. The present spectral interval does not contain suffi- 
cient spectral signatures of hot stars and thus cannot constrain 
the presence of hot stars in the solution. For the same reason, 
the presence of a featureless continuum in a composite spec- 
trum results in a high contribution of hot stars in the synthesis 
solution. 

Appendix A: The synthetic surface and the 
gradients of functions Q r : 

The "synthetic surface" is defined as the set of vectors of 
the normalised intensity vector space for which there exists a 
solution. In other words, the "synthetic surface" is the set of 
vectors I sy „ satisfying the equation system: 



Ilk - I sm = 
bk - 1 = 

The previous system can be written as follows : 

II ++ k + = 



(A.l) 



(A.2) 



where // + 



w 



and k + - 



1 syn 

i 

This equation is true only if for each combination of n* + 1 
lines of the matrix // ++ , the determinant of the resulting 
matrix equals zero (because vector k + is never equal to zero). 
Thus, as the rank of matrix 77 ++ is n+, one only needs to 
search for the n* linearly independent lines and impose that 
the determinants of the matrices obtained using these lines and 
to which is added a remaining line r (//,t + ), are equal to zero. 
Let us call the obtained determinants (which are of a number 
of «/ + 1 - n*) Q r where r is the index corresponding to the 
ni + 1 — n+ remaining lines. It follows: 



Q r = detail) = j]{-ir* +X+i I S ynidet(Ilt +Kn * +1) ) (A3) 



where //,. 



++ /(«« + !) 



is the extracted matrix from matrix + by 



eliminating the i' h line and the last column (i.e. the + I th col- 
umn, which is vector I syn ). 

In order to shift the origin of the vector space to I syn o, we per- 
form a translation of vector -/.,,„ o to the matrices II + r + and 
obtain: 



ftiw =(-ir* +1+ ^ f (//; 

di \isynp ~ otherwise. 



++;(«*+!) 



) ifj = h 



(A.4) 
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Unconstrained 


Dec. IMF 


Standard 


Star 


solution (%) 


mode (%) 


mode (%) 


07-BOV 











B3-4V 


2+1.5 


2± 1 


2± 1.5 


A1-3V 





0.06 ± 0.3 


0.2 ± 0.002 


F2V 











F8-9V 











G4V 





4±5 





G9-K0V 





0.6 ± 8 





K5V 


5 ± 19 


7.5 ± 19 


3± 19 


M2V 


0.3 ± 1 


1 ± 1 


0.1 ± 1 


rGOIV 











rG5IV 


4± 17 


8+13 


5 ±17 


rKOV 











rK3V 


10 ±8 


6±9 


9.5 ±8 


rMlV 











G0-4IH 


63 ±3 


56 ± 1 


64 ±3 


wG8HI 











G9III 











K4III 











M0.5IH 











M4III 


5±1 


5±1 


5± 1 


M5III 











rG9III 











rK3III 











rK3IIIbis 


7+12 


7± 10 


7± 12 


rK5III 


3± 13 


3 ±12 


3± 13 


GOIab 











K4Iab 











M2Ia 











rG2Iab 











rKOII 











rK3Iab 











D 1 


0.22 ± 0.017 


0.22 ± 0.017 


0.22 ±0.015 


E(B-V) 


0.30 ±0.10 


0.30 ±0.10 


0.31 ±0.1 


Dispersion (tr) 


60 ± 340 


60 ± 340 


60 ± 340 



Table 6. Results of the spectral synthesis of the globular cluster G 170. Each column displays the stellar contributions with their 
standard deviations for three solutions: the unconstrained solution and the solutions of the two modes of constraints described in 
the text. D 2 is the synthetic distance and E(B-V) the reddening. The dispersion is given in km/s. 



Appendix B: Facilities of the code: 

The Synthesis code is written in such a way that several entities 
can be comfortably changed or included. Here are listed some 
of them: 

- All components which may compose the galactic spec- 
trum consist of independent files containing the pixel val- 
ues of the corresponding spectra. Thus, the choice of the 
presupposed composition of the studied spectrum is com- 
pletely free. One may use all possible stellar, cluster and/or 
even galactic spectra. One may also use Blackbody spectra 
and/or any wavelength function. 

This facility can be used to test for example additive con- 
tinua in the case of AGNs. 

- The user of the code can choose wavelength intervals (from 
one pixel to the whole spectrum) to be excluded from the 
synthesis procedure. This is useful for example in the case 



of spectral regions showing emission lines or not well cali- 
brated. 

- Different weights can be included easily in the input of the 
program for each of the used pixel intensities of the spectra. 
This may be used to give more or less importance to some 
spectral regions in the synthesis procedure. In this case, the 
synthetic distance has the following form: 

d 2 =ruii' gair r synj ) 2 (B.D 

where /' , . and /' . are the weighed galactic and synthetic 

intensities (I' galj = I galj yfPj and I' synj = kjji y/Pj). 

- Dust spectra can be included easily in the database like the 
stellar spectra. They are also constructed with the required 
starting wavelength value, pixel step and pixel number us- 
ing a dedicated code. A subroutine of this code allows the 
user to decide the shape of the function describing the dust 
emission. 
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Unconstrained 


Dec. IMF 


Standard 


Star 


solution (%) 


mode (%) 


mode (%) 


07-BOV 











B3-4V 











A1-3V 











F2V 











F8-9V 











G4V 





16 ± 1.5 





G9-K0V 





3± 14 





K5V 





4± 11 





M2V 





1 ±0.4 





rGOIV 











rG5IV 











rKOV 


29 ±2 


3± 19 


29 ±2 


rK3V 


9±4 





9±4 


rMlV 











G0-4IH 


12+ 12 


7± 1.5 


12 ± 12 


wG8III 











G9III 





6± 12 





K4IH 











M0.5III 











M4III 


4± 1 


4 ±0.5 


4± 1 


M5III 


2 ±0.5 


2 ±0.5 


2 ±0.5 


rG9III 


45 ±6 


54 ±5 


45 ± 6 


rK3in 











rK3IIIbis 











rK5III 











GOIab 











K4Iab 











M2Ia 











rG2Iab 











rKOn 











rK3Iab 











D 1 


0.16 ±0.015 


0.17 ±0.015 


0.16 ±0.015 


E(B-V) 


0.09 ± 0.06 


0.08 ± 0.03 


0.09 ± 0.06 


Dispersion (tr) 


140 ± 60 


140 ± 60 


140 ± 60 



Table 7. Synthesis solutions for the globular cluster G 170 obtained using only the red part of the spectrum (from 6200 to 9000A). 



- The reddening law may also be changed easily, it consists 
of an independent file listing the pixel values of the law in 
the wavelength domain. 

- The dispersion is tested by convolving the spectra of the 
database with a Gaussian of chosen parameters. Any other 
function as the model using the Gauss-Hermite series can 
be easily included as a subroutine in the code. 

- At present, the constraints have the form of large inequali- 
ties between the different numbers (or equivalently the dif- 
ferent contributions to luminosity at the reference wave- 
length) of stars or objects of the database. Therefore, any 
linear constraints may be chosen in the code. 

- The interval values for all parameters can be changed con- 
fortably. 

An independent code "Proginterpol.f90" allowing users to 
homogenise the different spectra (the galactic or stellar spectra 
or the spectra of dust, reddening or any wavelength function) 
is also available. This code rebins the spectra to a single step 
and makes them all begin (first pixel) at the same wavelength. 
It allows one as well to obtain spectra having the same number 



of pixels. 



Appendix C: The output of the program: 

The code produces 4 files: 

- The "solution" which consists of a list of the contribu- 
tions (and their standard deviations) of each object of the 
database, the reddening and the dispersion with their er- 
rors. In this file are also listed all the input parameters of 
the synthesis. 

- The observed spectrum dereddened with the best obtained 
value. The file consists on an ASCII file listing the pixel 
values of the spectrum. 

- The synthetic constructed spectrum resulting from the com- 
bination of the contributions of the objects in the database. 

- An information file listing all used parameters and for each 
iteration the synthetic distance, the reddening, the disper- 
sion and the contributions of hot stars. This file is used to 
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plot figures showing the behaviour of the synthetic distance 
and hot stars with the different parameters as Figs. [2] and |4] 
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